function [IRF_VAR_boot_org,IRF_LP_boot_org,IRF_VAR_boot,IRF_LP_boot,IRF_VAR_boot_full,IRF_LP_boot_full]=VARLP_seqblockbootstrap(Y_var,nlags,IRF_hor,n_boot,responseV,normalizeV,scale,fixblock,blocksize)

[~,Y_var_lags,Y_var_ad]=lag_variable(Y_var,nlags);
[T_var,N_var]=size(Y_var_ad);
[M_var]=size(Y_var_lags,2);


e_store=nan(T_var*2,N_var,IRF_hor-1);
Phi_store=nan(M_var*2,N_var,IRF_hor-1);
rhs_varlp_h_store=nan(T_var,M_var,IRF_hor-1);
IRF_VAR_boot_full=nan(N_var^2,IRF_hor,n_boot);
IRF_LP_boot_full=nan(N_var^2,IRF_hor,n_boot);
IRF_VAR_boot=nan(n_boot,IRF_hor);
IRF_LP_boot=nan(n_boot,IRF_hor);
IRF_VAR_boot_org=nan(n_boot,IRF_hor);
IRF_LP_boot_org=nan(n_boot,IRF_hor);

for j=1:IRF_hor-1 
    
lhs_var_h=lag(Y_var_ad,j-1);
lhs_lp_h=lag(Y_var_ad,-j+1);
rhs_varlp_h=lag(Y_var_lags,j-1);

%stack lhs
lhs_varlp_h_ad=[lhs_var_h(j:end,:); lhs_lp_h(1:end-j+1,:)];
rhs_varlp_h_ad=[rhs_varlp_h(j:end,:)];
rhs_varlp_h_store(j:end,:,j)=rhs_varlp_h_ad;


%Estimate Phi and residuals
Phi_varlp=kron(eye(2),inv(rhs_varlp_h_ad'*rhs_varlp_h_ad)*rhs_varlp_h_ad')*lhs_varlp_h_ad;
e_resid=lhs_varlp_h_ad-(kron(eye(2),rhs_varlp_h_ad)*Phi_varlp);
e_var_resid=e_resid(1:size(e_resid,1)/2,:);
e_lp_resid=e_resid(size(e_resid,1)/2+1:end,:);

Phi_store(:,:,j)=Phi_varlp;
e_store(j:T_var,:,j)=e_var_resid;
e_store(T_var+j:end,:,j)=e_lp_resid;

end


%Block Bootstrap
for i_boot=1:n_boot
    for j=1:IRF_hor-1
    T_boot=size(rhs_varlp_h_store(j:end,:,1),1);
    if fixblock==0
    block_index=blockindex(T_boot,j-1);
    elseif fixblock==1
    block_index=blockindex(T_boot,blocksize-1);
    end

    block_index_full=[block_index; block_index+T_boot];
    rhs_varlp_call=rhs_varlp_h_store(j:end,:,j);
    e_store_call=[e_store(j:T_var,:,j); e_store(T_var+j:end,:,j)];
    Phi_store_call=Phi_store(:,:,j);
    bigeye= [eye(N_var) zeros(N_var,N_var*(nlags-1))];

  
    lhs_varlp_boot=kron(eye(2),rhs_varlp_call)*Phi_store_call+e_store_call(block_index_full,:);
    lhs_var_boot=lhs_varlp_boot(1:T_boot,:);
    lhs_lp_boot=lhs_varlp_boot(T_boot+1:end,:);


    Phi_boot=kron(eye(2),inv(rhs_varlp_call'*rhs_varlp_call)*rhs_varlp_call')*lhs_varlp_boot;
    Phi_var_boot=Phi_boot(1:M_var,:);
    Phi_lp_boot=Phi_boot(M_var+1:end,:);


    Phi_var_boot_int=Phi_var_boot(2:end,:);
    Phi_lp_boot_int=Phi_lp_boot(2:end,:);
    
    F_var_boot = [Phi_var_boot_int';eye(N_var*(nlags-1)) zeros(N_var*(nlags-1),N_var)];
    F_lp_boot  = [Phi_lp_boot_int';eye(N_var*(nlags-1)) zeros(N_var*(nlags-1),N_var)];


    if j==1
    %Estimate contemporaneous matrix
    e_resid_boot=lhs_varlp_boot-(kron(eye(2),rhs_varlp_call)*Phi_boot);
    e_var_resid_boot=e_resid_boot(1:T_boot,:);
    e_lp_resid_boot=e_resid_boot(T_boot+1:end,:);
    
    sigma_var_boot=e_var_resid_boot'*e_var_resid_boot/(T_boot-M_var);
    sigma_lp_boot=e_lp_resid_boot'*e_lp_resid_boot/(T_boot-M_var);

    if scale==0
     A0_var_boot=chol(sigma_var_boot,'lower');
     A0_lp_boot=chol(sigma_lp_boot,'lower');
    elseif scale==1
     A0_var_boot=chol(sigma_var_boot,'lower');
     A0_lp_boot=chol(sigma_lp_boot,'lower');
       for i=1:N_var
           A0_var_boot(:,i)=A0_var_boot(:,i)./A0_var_boot(i,i);
           A0_lp_boot(:,i)=A0_lp_boot(:,i)./A0_lp_boot(i,i);
       end  
    end

    %Estimate IRFs at h=0 and 1
    
    %h=0
    IRF_VAR_boot_full(:,j,i_boot)= reshape(A0_var_boot,[],1);  
    IRF_LP_boot_full(:,j,i_boot)= reshape(A0_lp_boot,[],1); 
    %h=1
    IRF_VAR_boot_full(:,j+1,i_boot)= reshape(bigeye*F_var_boot*bigeye'*A0_var_boot,[],1);
    IRF_LP_boot_full(:,j+1,i_boot)= reshape(bigeye*F_lp_boot*bigeye'*A0_lp_boot,[],1); 


    elseif j>=2  
    IRF_VAR_boot_full(:,j+1,i_boot)= reshape(bigeye*(F_var_boot^(j))*bigeye'*A0_var_boot,[],1);
    IRF_LP_boot_full(:,j+1,i_boot)= reshape(bigeye*F_lp_boot*bigeye'*A0_lp_boot,[],1); 
    end

    end
  
end

%Org and Scale bootstrap IRFs
for i_boot=1:n_boot
    IRF_VAR_boot(i_boot,:) = IRF_VAR_boot_full(responseV,:,i_boot) / IRF_VAR_boot_full(normalizeV,1,i_boot);  
    IRF_LP_boot(i_boot,:) = IRF_LP_boot_full(responseV,:,i_boot) / IRF_LP_boot_full(normalizeV,1,i_boot); 

    IRF_VAR_boot_org(i_boot,:) = IRF_VAR_boot_full(responseV,:,i_boot);
    IRF_LP_boot_org(i_boot,:) = IRF_LP_boot_full(responseV,:,i_boot);

end





end